{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# crystal_space_group calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The crystal_space_group calculation style characterizes a bulk atomic system (configuration) by determining its space group by evaluating symmetry elements of the box dimensions and atomic position.  This is useful for analyzing relaxed systems to determine if they have transformed to a different crystal structure.  An ideal unit cell based on the identified space group and the system's box dimensions is also generated.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Function slightly updated\n",
    "- 2020-09-22: Setup and parameter definition streamlined. Method and theory expanded.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "- [spglib](https://atztogo.github.io/spglib/python-spglib.html)\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "\n",
    "- The results are sensitive to the symmetryprecision parameter as it defines the tolerance for identifying which atomic positions and box dimensions are symmetrically equivalent.  A small symmetryprecision value may be able to differentiate between ideal and distorted crystals, but it will cause the calculation to fail if smaller than the variability in the associated system properties.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "The calculation relies on the spglib Python package, which itself is a wrapper around the spglib library.  The library analyzes an atomic configuration to determine symmetry elements within a precision tolerance for the atomic positions and the box dimensions.  It also contains a database of information related to the different space groups.\n",
    "\n",
    "More information can be found at the [spglib homepage](https://atztogo.github.io/spglib/).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)\n",
    "\n",
    "- [spglib](https://atztogo.github.io/spglib/python-spglib.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import datetime\n",
    "from copy import deepcopy\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://atztogo.github.io/spglib/python-spglib.html\n",
    "import spglib\n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'crystal_space_group'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is generated using the load parameters and symbols."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('prototype', 'A1--Cu--fcc', symbols='Ni', a=3.52)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Specify calculation-specific run parameters\n",
    "\n",
    "- __symmetryprecision__ is a precision tolerance used for the atomic positions and box dimensions for determining symmetry elements.  Default value is '0.01 angstrom'.\n",
    "\n",
    "- __primitivecell__ is a boolean flag indicating if the returned unit cell is to be primitive (True) or conventional (False).  Default value is False.\n",
    "\n",
    "- __idealcell__ is a boolean flag indicating if the box dimensions and atomic positions are to be idealized based on the space group (True) or averaged based on their actual values (False).  Default value is True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "symmetryprecision = uc.set_in_units(0.01, 'angstrom')\n",
    "primitivecell = True\n",
    "idealcell = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. crystal_space_group()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crystal_space_group(system, symprec=1e-5, to_primitive=False,\n",
    "                        no_idealize=False):\n",
    "    \"\"\"\n",
    "    Uses spglib to evaluate space group information for a given system.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    system : atomman.System\n",
    "        The system to analyze.\n",
    "    symprec : float\n",
    "        Absolute length tolerance to use in identifying symmetry of atomic\n",
    "        sites and system boundaries.\n",
    "    to_primitive : bool\n",
    "        Indicates if the returned unit cell is conventional (False) or\n",
    "        primitive (True). Default value is False.\n",
    "    no_idealize : bool\n",
    "        Indicates if the atom positions in the returned unit cell are averaged\n",
    "        (True) or idealized based on the structure (False).  Default value is\n",
    "        False.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Results dictionary containing space group information and an associated\n",
    "        unit cell system.\n",
    "    \"\"\"\n",
    "    # Identify the standardized unit cell representation\n",
    "    sym_data = spglib.get_symmetry_dataset(system.dump('spglib_cell'), symprec=symprec)\n",
    "    ucell = spglib.standardize_cell(system.dump('spglib_cell'),\n",
    "                                    to_primitive=to_primitive,\n",
    "                                    no_idealize=no_idealize, symprec=symprec)\n",
    "    \n",
    "    # Convert back to atomman systems and normalize\n",
    "    ucell = am.load('spglib_cell', ucell, symbols=system.symbols)\n",
    "    ucell.atoms.pos -= ucell.atoms.pos[0]\n",
    "    ucell = ucell.normalize()\n",
    "    \n",
    "    # Throw error if natoms > 2000\n",
    "    natoms = ucell.natoms\n",
    "    if natoms > 2000:\n",
    "        raise RuntimeError('too many positions')\n",
    "\n",
    "    # Average extra per-atom properties by mappings to primitive\n",
    "    for index in np.unique(sym_data['mapping_to_primitive']):\n",
    "        for key in system.atoms.prop():\n",
    "            if key in ['atype', 'pos']:\n",
    "                continue\n",
    "            value = system.atoms.view[key][sym_data['mapping_to_primitive'] == index].mean()\n",
    "            if key not in ucell.atoms.prop():\n",
    "                ucell.atoms.view[key] = np.zeros_like(value)\n",
    "            ucell.atoms.view[key][sym_data['std_mapping_to_primitive'] == index] = value\n",
    "    \n",
    "    # Get space group metadata\n",
    "    sym_data = spglib.get_symmetry_dataset(ucell.dump('spglib_cell'))\n",
    "    spg_type = spglib.get_spacegroup_type(sym_data['hall_number'])\n",
    "    \n",
    "    # Generate Pearson symbol\n",
    "    if spg_type['number'] <= 2:\n",
    "        crystalclass = 'a'\n",
    "    elif spg_type['number'] <= 15:\n",
    "        crystalclass = 'm'\n",
    "    elif spg_type['number'] <= 74:\n",
    "        crystalclass = 'o'\n",
    "    elif spg_type['number'] <= 142:\n",
    "        crystalclass = 't'\n",
    "    elif spg_type['number'] <= 194:\n",
    "        crystalclass = 'h'\n",
    "    else:\n",
    "        crystalclass = 'c'\n",
    "    \n",
    "    latticetype = spg_type['international'][0]\n",
    "    if latticetype in ['A', 'B']:\n",
    "        latticetype = 'C'\n",
    "    \n",
    "    pearson = crystalclass + latticetype + str(natoms)\n",
    "    \n",
    "    # Generate Wyckoff fingerprint\n",
    "    fingerprint_dict = {} \n",
    "    usites, uindices = np.unique(sym_data['equivalent_atoms'], return_index=True)\n",
    "    for usite, uindex in zip(usites, uindices):\n",
    "        atype = ucell.atoms.atype[uindex]\n",
    "        wykoff = sym_data['wyckoffs'][uindex]\n",
    "        if atype not in fingerprint_dict:\n",
    "            fingerprint_dict[atype] = [wykoff]\n",
    "        else:\n",
    "            fingerprint_dict[atype].append(wykoff)\n",
    "    fingerprint = []\n",
    "    for atype in sorted(fingerprint_dict.keys()):\n",
    "        fingerprint.append(''.join(sorted(fingerprint_dict[atype])))\n",
    "    fingerprint = ' '.join(fingerprint)\n",
    "\n",
    "    # Return results\n",
    "    results_dict = spg_type\n",
    "    results_dict['ucell'] = ucell\n",
    "    results_dict['hall_number'] = sym_data['hall_number']\n",
    "    results_dict['wyckoffs'] = sym_data['wyckoffs']\n",
    "    results_dict['equivalent_atoms'] = sym_data['equivalent_atoms']\n",
    "    results_dict['pearson'] = pearson\n",
    "    results_dict['wyckoff_fingerprint'] = fingerprint\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = crystal_space_group(ucell,\n",
    "                                   symprec=symmetryprecision,\n",
    "                                   to_primitive=primitivecell,\n",
    "                                   no_idealize=not idealcell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Display space group information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number\n",
      "225\n",
      "\n",
      "international_short\n",
      "Fm-3m\n",
      "\n",
      "international_full\n",
      "F 4/m -3 2/m\n",
      "\n",
      "international\n",
      "F m -3 m\n",
      "\n",
      "schoenflies\n",
      "Oh^5\n",
      "\n",
      "hall_symbol\n",
      "-F 4 2 3\n",
      "\n",
      "choice\n",
      "\n",
      "\n",
      "pointgroup_schoenflies\n",
      "m-3m\n",
      "\n",
      "pointgroup_international\n",
      "Oh\n",
      "\n",
      "arithmetic_crystal_class_number\n",
      "72\n",
      "\n",
      "arithmetic_crystal_class_symbol\n",
      "m-3mF\n",
      "\n",
      "ucell\n",
      "avect =  [ 2.489,  0.000,  0.000]\n",
      "bvect =  [ 1.245,  2.156,  0.000]\n",
      "cvect =  [ 1.245,  0.719,  2.032]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 1\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "\n",
      "hall_number\n",
      "523\n",
      "\n",
      "wyckoffs\n",
      "['a']\n",
      "\n",
      "equivalent_atoms\n",
      "[0]\n",
      "\n",
      "pearson\n",
      "cF1\n",
      "\n",
      "wyckoff_fingerprint\n",
      "a\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for key in results_dict.keys():\n",
    "    print(key)\n",
    "    print(results_dict[key])\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
